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We use a particle-based mesoscale model that incorporates chemical reactions at a 
coarse-grained level to study the response of materials that undergo volume-reducing 
chemical reactions under shockwave-loading conditions. We find that such chemi¬ 
cal reactions can attenuate the shockwave and characterize how the parameters of 
the chemical model affect this behavior. The simulations show that the magnitude 
of the volume collapse and velocity at which the chemistry propagates are critical 
to weaken the shock, whereas the energetics in the reactions play only a minor role. 
Shock loading results in transient states where the material is away from local equilib¬ 
rium and, interestingly, chemical reactions can nucleate under such non-equilibrium 
states. Thus, the timescales for equilibration between the various degrees of freedom 
in the material affect the shock-induced chemistry and its ability to attenuate the 
propagating shock. 


I. INTRODUCTION 

When materials are subject to dynamical mechanical loads (shockwaves) a 
plethora of complex processes are launched in response to the insult. These can 
include plastic deformation, Una, phase transitions, M, chemical reactions mini, 
and even electronic transitions |^. Studying the response of different materials to 
shockwaves has resulted in significant insight into the response of materials under 
extreme conditions. na In most cases, when solid materials are shocked above a 
threshold strength, known as the Hugoniot elastic limit, plastic deformation nucle¬ 
ates behind the shock front to release the compressive stress along the shock direction 
and minimize free energy by achieving a more hydrostatic state. This stress 

relaxation weakens the leading shockwave, often resulting in a two-wave structure, 
with a leading elastic wave followed by a plastic wave that propagates at slower 
speeds; similar two-wave structures have been found in shock induced martensitic 
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transformations |1]. In the case of explosives, exothermic reactions that lead to low 
density products can enhance the initial shockwave and transform it into a deto¬ 
nation wave. For this to happen, the chemical reactions have to generate enough 
energy and pressure. However, thermodynamics alone is not sufficient; fast reaction 
kinetics are critical to significantly affect a traveling shockwave. 

Materials that dissipate or absorb the energy in shockwaves causing the shock front 
to weaken as it propagates are of interest for protection against impacts, collisions and 
blasts. Impact resistant materials include fiber composites |I3] and ceramic/metal 
armor [H]. There has been growing interest in softer materials capable of absorbing 
shockwave energy via molecular-level processes. Polyurea, a polymer with glass 
transition below room temperature, has been shown effective at absorbing energy 
in ballistic impact tests. ng While the mechanisms for the high energy absorption 
capability at high strain rates are not fully understood, a transition to the glassy 
state is believed to play a key role. dSHni Materials with very high porosity at the 
nanoscale, such as metal organic frameworks, are also attracting attention for such 
applications as void collapse can contribute to weaken the shockwave [18] . 

In this paper, we explore the possibility of using chemical reactions in molecular 
systems for shock wave energy absorption; specifically endothermic, volume reducing 
chemistry. We use a recently developed model for coarse grain simulations for a 
class of model materials that exhibit the desired behavior at the molecular level 
and and study the mechanisms through which shock energy dissipation can occur. 
Our objective is to understand requirements of such shockwave energy dissipating 
(SWED) materials and to establish how the characteristics of the chemical reactions 
affect its ability to weaken shocks. We foresee this knowledge could contribute to 
experimental design efforts. Our simulations demonstrate that a chemical reaction 
front involving endothermic, volume-reducing chemistry can propagate fast enough to 
couple with the leading shockwave and weaken it. The results also shed light on how 
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the nucleation and propagation of such chemical reactions occur under dynamical 
loading. 

The paper is organized as follows. Section II introduces ChemDID, the coarse 
grain model used to investigate shock-induced chemistry and the family of model 
materials to be characterized. Section III discusses the response of materials that 
can undergo volume reducing chemistry to shock loading and relate characteristics 
of the chemical reactions and the ability of the materials to weaken shocks. Section 
IV discusses the molecular processes that govern the nucleation of chemistry behind 
the shockwave and Section V draws conclusions. 


II. MODEL AND SIMULATION DETAILS 

A. ChemDID 

ChemDID [Tn| is a coarse-grained, particle-based model that enables the descrip¬ 
tion of stress-induced chemical reactions involving degrees of freedom (DoF) internal 
to the mesoparticles. Unlike all-atom MD, where particles represent atoms, in Chem¬ 
DID particles describe extended objects, molecules in this paper. Such molecules will 
be represented by a spheres with 3N-3 internal DoFs (for an N atom molecule), out 
of which we will single out one of them to describe the chemical reactions in the 
molecule while the rest will be treated via the equipartition theorem. 

The variable singled out will represents the molecular radius (a), and is described 
with explicit Hamiltonian dynamics and an associated potential energy that can 
enable chemical reactions. The remaining 3N-4 DoFs are described in an averaged 
form using the approach proposed in Ref. ESI. Their state is described by a single 
dynamical variable that represent their temperature. Thus, the complexity of many- 
body intra-molecular interactions, including the desired volume collapsing reactions. 
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Intra-molecular 

Potential 


FIG. 1. Depiction of intermolecular and intramolecular potentials in ChemDID 

is reduced to an Intra-molecular potential inspired in transition state theory, where 
the size of the mesoparticle is used as an order parameter that governs the transition 
between an initial high-volume low-energy state or a collapsed-low volume high- 
energy state; we will refer to such materials as SWED materials, while materials 
that do not undergo chemical reactions will be called Inert. Fig. [^contrasts the two 
cases described. The non-bonded interactions between mesoparticles are described 
by pair-potentials which acts from surface to surface distance as shown in Figure 

The Hamiltonian of the system is: 


{cTj}, {Pi}, {tT*}) = y; (pinter{\Ri “ Rj 


- Gi- a 



i<j 



( 1 ) 


where (f>inter and (f>intra describe the inter-molecular and intra-molecular potential 
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radius (cr) [A] 


(a) (b) 

FIG. 2. (a) Representation of ChemDID as a coupled spring system with an intermolecular 
and intramolecular terms, (b) The Intra-molecular potential used describes an Inert case 
(dashed line) and SWED case (solid line); The parameters AG determine the activation 
barrier energy and AH describes the amount of energy absorbed during a collapse to a 
low-volume state. 


terms, Pi describes the translational momentum and associate mass m of particle i, 
and TTi similarly describes the conjugate momentum to the breathing mode with its 
associated inertial parameter m*. A derivation of the equations of motions has been 
previously shown in ra 

We note that the Hamiltonian does not account for the remaining 3N-4 internal 
DoFs of the molecules. These modes exchange energy with the Hamiltonian vari¬ 
ables (center of mass position and molecular radius) and they are described with the 
approach proposed in Ref. [20]. These internal DoFs are incorporated statistically 
and described by temperature (T*”*); they are coupled to the explicit DoFs via the 
position update equation such that energy flows to equalize the temperatures asso¬ 
ciated with the various degree of freedom: that of the radial breathing mode 

Trad q£ q£ molecules . The resulting equation of motion 
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for the internal temperature, see Refs. [201 EU , is: 


rpint 

_ rjiint _ , 

(Jint ~ ~ ^rneso 


(tt 


j^int ^ 






Hnter \ 2 


^rad 


^rprad _ rpint^ 

m\Cf\ulJQc 


prad 12 


( 2 ) 


where is the energy of the implicit DoFs, Q”* is their specific heat, r'meso 
and Urad describe the strength of the internal-to-intermolecular coupling and the 
internal-to-radial coupling respectively, 0o is a reference temperature, and the ratio 
|.Fp/(m(a;^)) provides a natural timescale for the corresponding interaction. Note 
that the magnitude and direction of the energy flow between internal DoFs and the 
molecular centers of mass (first term in the RHS of Eq. and the breathing mode 
(second term) is governed by the difference in local temperatures such that heat flows 
from hot to cold. The total energy is a conserved quantity in ChemDID and it is 
given by 


Etot ^ ^ Winter (IRi Rj I T ^ ^ 4^intra(^(^i) 



^k^T^meso 


( 3 ) 


We remind the reader that although the above definition of temperatures only ap¬ 
plies for equilibrium conditions after time or ensemble averages in the canonical 
ensemble, defining instantaneous and local values is useful to study processes out-of- 
equilibrium. 
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B. A model reactive molecular crystal 


The parameterization of a ChemDID model involves determining inter-molecular 
and intra-molecular potentials, the molecular mass and the inertial parameter for 
the dynamics of the molecule radii and as well as the coupling constants, I'meso and 
I'rad, that describe the coupling between the internal DoFs and the molecular centers 
of mass and radii. In this paper we parameterized a model reactive material with 
thermo-mechanical properties similar to anthracene, a molecular material believed to 
be capable of endothermic, volume reducing chemistry; DFT calculations by Slepetz 
et al. [22] show that the low-volume endothermic state of anthracene has an energy 
between 10-20 kcal/mol over the high-volume ground state. 

We studied a family of intramolecular potential associated with the breathing 
mode in order to quantify how characteristics of the chemical reactions affect the 
shock weakening power of the material. We characterize the chemical reactions 
by three key parameters: i) volume collapse, ii) endothermicity, and hi) activation 
energy. The initial radius of the molecule is taken as cr = 2.50 A, this value corre¬ 
sponds roughly to one of the dimensions that undergoes the greatest change in the 
anthracene molecule under the application of pressure [221 EHj- The intra-molecular 
potential used takes the form: 


(t>intra{<y) = K * {a - Umin) ' {<7 “ <7max) + Ai? - - _ . + Ae 

[^min ^max) 

( 4 ) 

,where (Jmin and a^ax denote the low-volume and high-volume (meta)stable points, 
the parameter K determines the overall curvature of the potential away from the 
(meta)stable points, while a Gaussian term (A,do) controls the curvature of the 
region in-between, /S.H denotes the amount of endothermicity of the reaction for the 
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collapsed state at cr = cTmm, and the activation barrier AG will depend implicitly on 
both the Gaussian parameters (A,cro) and the overall curvature parameter K. Table 
H] shows the intramolecular potential parameters. 

In determining the overall collapse of the molecule, the radius of the sphere and 
the inter-molecular separation will both play a role. In this paper, we quantify this 
volume-collapse in terms of its van der Waals radius 

4 

Vvdw = ^7r( (T -|- Aydw) ^ (5) 

f'vdW 

where van der Waals radius (r^dvi/) will be given by the sum of the sphere’s hard-core 
radius (cr) and a van der Waals skin {Ay^w)- The later depends on the intermolecular 
potential, discussed next, and takes the value to 1.07A. A Morse potential will be 
used to describe its inter-molecular interactions, given by: 

(t>inter{Rij - (Ji - ^j) = .-d/ro)] (6) 

where the parameters cq, tq, and 7 describe the cohesive energy, interaction range, 
and its curvature near its global minima; the values used are shown in Table |I} The 
ground state structure of such a system is the fee crystal structure and the inter- and 
intra-molecular parameters chosen resulting a lattice parameter of a = 10.1 A. The 
molecular mass will be taken asm = 296.1 g/mol, which has been used previously as 
a benchmark for molecular crystals and corresponds to an HMX-molecule [121 El] • 
This gives an equilibrium volume per molecule (Wg ) of 258.8 A^, close to that of 
athracene (231.2 A^). The equilibrium volume per molecule and the van der Waals 
volumes are related by the packing fraction (= which corresponds to a value 

of 0.74 for an fee arrangement. 

The remaining parameters in the model determine the coupling constants between 
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Vibrational Frequencies vs Inertial Parameter 



(a) (b) 


FIG. 3. (a) Black and red curves correspond to the phonon and radial breathing modes 
- with different m*- respectively vs frequency, (b) Dotted and solid lines correspond to 
the mean vibrational frequency squared (cj^) for the phonon and breathing modes (for 
different volume-collapse) respectively vs inertial parameter (m*); the expected behavior 
(cj^) oc ^ can be observed, where the proportionality factor is given by the curvature of 
the intramolecular potential at its global minima. 


the various DoFs and to do this it is useful to look at ChemDID the in terms of a 
coupled spring model, where a set of intra-molecular spring are in series with inter- 
molecular springs as shown in Fig. [^(a). Similar vibrational frequencies between the 
intra-molecular and inter-molecular modes allow for an effective exchange of energy 
between them, while a significant difference will tend to decouple the two sets of 
modes. For this reason the rate of energy exchange between the translational mode 
l^jimeso^ and radial modes (T^^^) can be tuned by changing the inertia parameter(m*) 
which increases the overlap between the intra-molecular (sphere modes) and the 
translational phonon spectrum, and therefore, leads to a stronger coupling. Figure 
(a) shows the (un-normalized) vibrational spectrum for the breathing and the phonon 
modes for different values of the inertia parameter (m*). Figure]^ (b) shows the mean 
squared frequency (cj^) = f {27TuYP{u)dh'^ where P{iy) is the normalized density of 
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states in Figure]^ (a). Note that the value of m* is chosen as to bring the mean 
of the sphere vibrational frequency to the same range as the mean of the phonon 
vibrational frequency. Table [I] shows how these values depend on the compliance 
between the two springs. 


III. SHOCK PROPAGATION AND ATTENUATION 

A. Response to a sustained shock 

In the following we compare the shock response of inert and SWED materials 
using ChemDID. The initial condition for the simulations consist of a target made 
of an fee crystal obtained by replicating the four-atom fee unit cell (with lattice 
parameter of 10.1 A) 200 times along the shock direction (z) and 20 times along the 
X and y directions leading to a system with 320,000 molecules. Periodic boundary 
conditions are imposed along the x and y directions, while the z direction is open. 
The system is thermalized at 300 K for 100 ps bringing internal DoFs in ChemDID 
to thermal equilibrium with radial and center of mass modes. 

The samples are impacted with a thin (one lattice constant thick), rigid and 
infinitely massive piston traveling at the piston velocity (up). For piston’s speeds 
below 2.0 km/sec we use an integration time step of dt = 0.005 ps, whereas for 
speeds above this value we halve the time step in order to get numerical stability. The 
infinitely massive piston does not slow down due to the interactions with the target, 
nor rarefaction waves are generated from its free surface. Thus this setup generates 
a sustained shock as opposed to one of finite duration obtained with finite pistons. 
Since the piston does not slow down, the material that is in its immediate proximity 
will quickly approach the same velocity Up. The material under stress will attempt 
to deform or rearrange in order to alleviate its local stress build up. While in an inert 
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TABLE 1. Parameters used in ChemDID 

Internal Parameters 

parameter symbol unit value 


Internal Heat 


Capacitance 

Cint 


— 

60 

internal-to-intermolecular 





coupling 

^meso 

l/ps 

0.1 

internal-to-intramolecular 





coupling 

^intra 

1/ps 

0.1 

mesoparticle mass 

m 


g/mol 

296.16 

Intermolecular parameters 

parameter 

symbol 

units 

values 

morse potential 

Winter (-Pij 


, 7 ,ro) kcal/mol 

Eqn. lil 

range 

To 


A 

2.14 

curvature 

7 


— 

4.5 

energy 



kcal/mol 

7.0 

mean freq. square 

inter) 

l/ps^ 

120.0 

Intramolecular parameters" 

1 

parameter 

symbol 

unit 

values 

ChemDID potential 

Sintra i ^min : 

AG, AH) kcal / mol 

Eqn. 1^ 

barrier 

AG 


kcal/mol 

30 - 80 

endothermicity 

AH 


kcal/mol 

0 - 20 


parameters 
depending on amin 

^min 

1.50 A 

1.75 A 

2.00 A 

2.25 A 

volume change 

(AV/VoUw 

— 

65 % 

50 % 

35 % 

20 % 

inertial parameter 

m* 

g/mol 

2961 

4145 

7402 

21319 

mean freq. 


l/ps2 

118.4 

123.00 

120.3 

119.2 

curvature 

K 

— 

80 

253 

1280 

20500 

Gaussian width 

CTo 

A 

0.177 

0.133 

0.088 

.044 


Gaussian constant A kcal/mol AG — — 5.0 

maximum radius cFmax ^ 2.5 


^ Parameters here have been defined such that AG > AH and AG > 5 [kcal/mol] 
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FIG. 4. Snapshots showing (vdW) volume, temperature, and crystal structure for INERT 
(left) and SWED (right) materials impacted from the left by an infinitely massive piston 
with velocity of 2 km/sec 


material, plastic deformation or phase transformations are the only mechanisms that 
allow for stress relaxation, in a SWED material volume-reducing chemical reactions 
can significantly reduce pressure build up and attenuate the leading shock, as we 
shall see. 

Figure]^ shows atomic snapshots of steady shocks on a SWED and inert materials 
shocked with Up=2.0 km/s. Molecules are colored to indicate local volume, temper¬ 
ature, and crystal structure. Interestingly, the SWED material exhibits a region of 
reacted (volume collapsed) material next to the piston; this reacted region propa¬ 
gates along the shock direction. As expected, in the inert case shocked particles 
also reduce their volume due to the high pressures but this occurs rather uniformly 
throughout the sample. Importantly, we see that the leading shock has propagated 
much further in the inert material (all snapshots correspond to the same time of 25 
ps), an indication that the chemistry is weakening the shock. The amount of plastic 
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deformation following the shock also sheds light into its strength. We analyze plas¬ 
tic and structural transformations by characterizing the local crystal environment of 
each molecule. Hep atoms in the fee crystal indicate stacking faults that separate 
partial dislocations; thus red atoms indicate the traces of partial dislocations or re¬ 
gions that transformed to hep. We see significantly less plastic deformation in the 
SWED material; this is because the leading shock is weaker. A temperature spike 
follows the primary wave in the inert case, but it in the SWED case the highest 
temperatures localize in the chemically reacted region. 
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FIG. 5. Velocity, Temperature, Pressure , and density profiles for a inert (solid) and swed 
(dashed/dotted) samples. The piston speed is Up = 1.25 [km/sec] 


Figure [^contrasts the profiles of key thermodynamic quantities for a SWED and 
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FIG. 6. (LEFT) Velocity profiles at different times, (RIGHT) Position of the shock fronts 
vs time 

inert materials at a given time; the SWED profiles further contrast the response 
with zero {AH = 0 kcal/mol) and finite endothermicity {AH = 20 kcal/mol). A 
three-wave structure can be seen for the SWED material. The leading wave (often 
called elastic precursor) shocks the material and behind it a plastic wave propagates 
at a slower velocity. Following the plastic is a chemical wave that densifies and heats 
up the material. Note that higher endothermicity decreases the temperature in the 
reacted region. Also the stiffness here will be slightly higher due to the linear term 
in AH appearing in the intramolecular potential. It is interesting to note that the 
pressure is greater in the reacted material than in the rest of the sample. It will be 
shown below, that this is a consequence of the volume collapse during a sustained 
shock and a necessity for steady state shocks. A similar three-wave structure has 
been observed experimentally on porous copper |25] . 

Figure]^ (a) shows velocity profiles at different times and depicts the development 
of the three wave structure. The chemical wave separates from the initial shock 
at short times and the plastic and elastic waves take significantly longer time to 
separate. For the case shown the separation between the elastic and plastic waves 
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occurs at around 30 ps, after which the two waves trail each other within a few 
nm. In contrast, the inert material is in an overdriven regime where the plasticity 
has merged with its elastic propagation into a single wave. Figure]^ (b) shows the 
position of the three wave fronts as a function of time from which we obtain their 
velocities. We see that steady-state is reached very early for the chemical wave, 
whereas for the plastic/elastic fronts it takes a longer time for this to occur. 


B. Taming shockwaves with volume-reducing chemical reactions 

Now that we have established the possibility of weakening shocks with chemical 
reactions that involve volume collapse, we are interested in understanding how the 
characteristics of the chemistry affect the strength of the coupling and the SWED 
effect. Before describing our ChemDID results we use conservation laws to develop a 
general framework to discuss SWED and understand how volume collapse, activation 
barrier and endothermicity contribute to achieving the desired effect. 

We start by considering the locus of the pressure-volume states accessible by 
shocking a material initially at volume Vq and negligible pressure Pq ~ 0, see blue 
dots in Figure It is important to understand that these states (the Hugoniot 
of the material) is different from a reversible equation of state; each point is the 
result of a shock experiment and can only be accessed by shocking material at the 
initial point (Vo,.Po)- Assuming the shock is too weak to trigger chemistry or plastic 
deformations, we will have a single wave structure m- By applying mass and 
momentum conservation across the shock front we will obtain interesting insight 
into the SWED behavior. 

Requiring mass conservation across the shock wave [2S] yields the following equal- 
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FIG. 7. Schematic showing the paths in the P-V space for inert and reactive shocks. The 
critical Hugoniot Chemical Limit (HCL) connects the transition point between the two 
Hugoniot curves. 


ity: 


ps{Xs - Us) = - Us) (7) 

where pa and Xg represent the density and particle velocity in the shocked region and, 
as before, snbscript 0 denote the unshocked region. Since the unshocked material is 
not moving {xo = 0), we can readily solve for the particle velocity in the shocked 
region in terms of the densities and the shock velocity: 

= ( 8 ) 

Ps 


An expression for the pressure in the shocked region can be written from the 
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momentum conservation equation as: 


Ps = poUsXs = Poulil -—) (9) 

Ps 

The case when the pre-shocked region is not stationary is given in the appendix 
VI A[ Rearranging the two conservation equations, known as Rankine-Hugoniot 
equations, m we can write down an expression for the shock velocity in terms of the 
change in pressure and volume 


2 2 


Ps-Po 
Vs-Vo 


( 10 ) 


We see that the slope of the line connecting the initial and shocked materials 
(black lines in Fig. known as Rayleigh lines) is related to the velocity squared of 
the corresponding wave. It is clear from the conservation equations that in order 
to decrease the pressure following the initial shock it is necessary to slow down the 
shock speed Ug or, equivalently, maintain the density in the shocked region as close 
as possible to the un-shocked density. Let’s delve into how chemical reactions can 
help to maintain both of these conditions. 

If we shock the material above a given threshold, that we will call Hugoniot 
chemical limit (HCL) in analogy to the Hugoniot elastic limit, a chemical wave will 
develop behind the elastic one. For simplicity we will consider two wave structures 
and neglect the plastic wave. In the pressure-volume plot we connect the state 
representing the elastic wave with the Hugoniot of the chemical products (red circles 
in Fig. [^. We can see that there will be a series of shocks for which the velocity of the 
chemical wave will be lower (gentler slopes) than that of the leading elastic precursor 
(steeper slope). In this regime, the point {Phcl,Vhcl) in the Hugoniot determines 
the propagation velocity and pressure of the elastic precursor through Eqn. (j^ and 
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(10), where s = HCL (assuming no plasticity). This value is independent of the 
piston velocity (the dependence of the HCL with the activation barrier is described 
in the next section) until the velocity of the chemical wave matches that of the elastic 
precursor. At this point, the overdriven regime is reached, and the prohle of the wave 
is characterized by a single wave structure. 


Let us now quantify how the chemical wave weakens the leading shock wave. 
Applying mass conservation across the chemical wave { c }, i.e. between the reacted 
and shocked regions leads to: 


Pc(^C ^c) Psi^s ^c) 


( 11 ) 


where pc and Xc represent the density and particle velocity in the reacted region. 
Substituting Eq. into Eq. and rearranging terms we obtain an expression for 
the shock speed as: 


Xg p i^c ^c) T 

7/ = - = — - 

^ I — Po I _ Po 

Ps Ps 

Moreover, in a sustained shock the particle velocity of the reacted material is equal 
to the piston velocity {xc = Up). Hence, we arrive at an expression for the pressure 
following the leading shock in terms of the densities of the various regions and the 
velocity of the chemical wave: 



Pg Po^S^S Po 


'^c) + '^c) 


1 


Po_ 

ps 


(13) 


Let us comment on the regime under which this expression is valid. In order 
for waves to separate into a chemical and a shocked components, the velocity of 
the chemical wave needs to lie between the piston velocity Up and the wave speed 
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corresponding to the HCL: . The upper limit, Uc = corresponds to the 

overdriven regime when chemical wave and the shock wave travel at a single speed 
n* = Yzfe- The lower limit, i.e. a chemical wave traveling at Uc < Up corresponds to 

Pc 

no chemistry production, viz. Uc —t 0 ^ pc —t therefore the (inert) shocked wave 

travels at speed: uf = 

Ps 



(a) (b) 

FIG. 8. (a) Expression [T^ (for Up=1.75 km/sec) compared to simulated data with differ¬ 
ent volume and chemical speeds; for a given volume collapse, a lower activation barrier 
corresponds to a faster chemical speed. The point where the chemical wave reaches the 
overdriven regime is denoted by (*). (b) Shock speed vs impact speed for various volumes 
(for the same activation barrier AG = 30 kcal/mol); The chemical waves are denoted by 
the dashed lines. 

Equationshows that in order to weaken the shock, i.e. reduce Pg, we need pc to 
be high (large volume collapse) but also we need the chemical wave to travel at fast 
speeds. This makes it clear that, as discussed earlier, reaction kinetics are critically 
important for SWED. We now compare the predictions of Eq. [T^ for the pressure 
behind the leading shockwave with results from explicit ChemDID for a family of 
SWED materials with various volume collapse amounts and activation energies. The 
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density in the shocked regions ps depends slightly on the activation barrier, but we 
will assume it to be constant at pg = 2500 kg/m^ to evaluate Eq. 13; this number 
changes at most by 20 % for all the cases considered at this piston speed. We shall 
see that this assumption leads to reasonable predictions and, in the next section, we 
will explore the sensitivity of the results with respect to it. 

Figure (a) shows the pressure following the initial shock as a function of the 


velocity of the chemical wave. The lines show the predictions from Eq. and the 
symbols represent ChemDID results. The results are shown for four model materials 
with different amounts of volume collapse and various activation barriers (AG = 
{80, 60,30} kcal/mol). As expected, increasing the amount of volume collapse results 
in lower shock pressures. The non-linear dependence of pressure on the density of 
the compressed state is consistent with the analysis of shocks on porous materials 
[2ZII2H]. For a given volume collapse, lowering the activation energy (that controls 
reaction kinetics) leads to faster chemical waves and also result in a reduction of 
the shock pressure. Figure (b) shows the effect of volume-collapse on the leading 
shockwave (open symbols) and chemical wave (filled symbols) speed as a function of 
the piston velocity; we show the same volume-collapse cases in (a). Higher volume 
collapse result in lower shock velocities for all piston velocity and a reduction of 
the velocity corresponding to the HCL. Once a chemical wave propagates, it couples 
with the leading shock wave resulting in a plateau of the shock wave speed as a 
function of the piston speed. Interestingly, we find that chemical waves tend to 
propagate faster the smaller the volume-collapse for a given piston velocity. These 
results are also consistent with the shock experiments on porous copper of different 
densities|2S]- The authors report a three-wave structure, similar to what is observed 
in our simulations, where measurement of the slower wave (equivalent to our chemical 
wave) shows a steeper slope for samples with less porosity. The point where the 


chemical wave merges with the leading shock wave can also be seen in Figures 
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K a) and (b) and is denoted by the asterisks. This points marks the end of the shock 
absorbing capabilities of the material and further increasing the piston velocity leads 
to increased shock velocity and pressure. 

In summary, we have shown that volume collapse is the dominant factor in deter¬ 
mining the effectiveness of the chemical reactions at weakening the leading shockwave 
(i.e. reducing the pressure) but reaction kinetics is also important. On the other 
hand, large volume collapses cases exhibit HCLs corresponding to slower shock ve¬ 
locities and the design of materials for shock absorption should take into account the 
regime of operation. 


C. Quantifying the role of chemistry on shock propagation 

As discussed in the previous sub-section, in order to engineer materials with de¬ 
sired shock-wave dissipation response, it is imperative to establish relationships be¬ 
tween the characteristics of the chemical reactions and the material response. We 
have already shown the importance of volume collapse and the velocity of the chemi¬ 
cal wave. In this section we explore ChemDID simulations in more detail, in order to 
correlate microscopic characteristics of the chemical reactions with shock dissipation. 

To understand the development of a two-wave structure and its properties we 
show the evolution of the system in the pressure-volume space. Using ChemDID, 
we follow the local properties of a thin slab of material during shock loading. The 
results are shown in Fig. (a) for three different activation barriers for a volume 
collapse of (AU/Vo)„dvF = 35 %. This hgure shows the individual paths for a range 
of piston speeds between Up = 0.25 — 3.0 km/sec (in steps of 0.25 km/sec). For 
velocities below a critical value (HCL-limit; see Fig [^, the ends of the Rayleigh 
lines follow the inert (unreacted) Hugoniot. After a critical compression is reached, 
corresponding the volume of the HCL (= 1/p^^^) which depends slightly 
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(a) (b) 

FIG. 9. (a) Hugoniot and (b) Us-Up curves for an Inert and SWED materials with a 
curvature on the intramolecular potential corresponding to {ISV/VQ)ydw = 35 %. The 
different symbols correspond to different barriers. A single impact speed Up= 1.75 km/sec 
has been emphasized in (a) with larger symbols and thicker lines. The chemical waves are 
shown in dashed lines in (b). 


on the activation barrier for chemistry, the system develops a two-wave structure 
described by two Rayleigh lines: the first going from the initial state to the pressure 
and volume after the initial wave (along the unreacted Hugoniot) and the second 
crossing into the reacted Hugoniot. In this description we are combing the plastic 
and elastic waves into one as their propagation velocities are similar and they do not 
separate significantly during the simulation time. 


In order to investigate the dominant factors during this transition, we can write 
down the following expression relating the pressures in the reacted region Pc and 


unreacted shocked region the derivation is shown in Section VIA of the Ap- 
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pendix. 
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The pressure in the reacted zone Pc will lie in the reacted Hugoniot and will depend 
primarily on the impact speed. Higher pistons will cause the reacted region to move 
towards higher pressures along the product Hugoniot. This can be seen in Figj^where 
the case Up=1.75 km/sec has been emphasized. The pressure response of the reacted 
material for this case will tend to collapse to the same region around y/Vo = .42 and P 
= 17 GPa. The critical pressure Pf^^, on the other hand, will depend on the onset of 


chemistry. As seen from Eqn. 14, this pressure is mainly dominated by density of the 
reacted material and the propagation velocity of the chemical wave since they appear 
quadratic in this expression, while the critical volume appears linearly. This 


is the reason we were able to ignore the crucial density in Eqn. 13 and obtain 

and expression that matches well with the simulated data. We also see from Eqn. 


14, that the change in pressure cannot be negative, since pc > pf"^- The slope 


of the Raleigh line in the P-V plane is therefore zero near threshold, and negative 
for case with more volume collapse. The slope of the Raleigh-line determines the 
propagation velocity of the chemical front. The points that lie between the reacted 
and inert Hugoniots represents cases with heterogeneous regions of collapsed volume 
which do fully propagate, at least within the scales of our simulation. 


In Eig. 1^ (b) we compare the shock velocities for an inert an reactive ChemDID 
simulations with various activation energies and includes experimental data on an¬ 
thracene for reference [22] • In the elastic regime, i.e. shocks weaker than the HEL, 
the shock speed increases linearly with piston velocity. The HEL is reached for piston 
velocities of approximately 0.75 km/sec, this value depends weakly on the stiffness 
of the intra-molecular potential; stiffer materials tend to have slightly higher HEL 
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limits. The HCL, on the other hand, occurs at slightly higher piston velocities and 
depends strongly on the activation barrier of the intra-molecular chemistry (dashed 
lines) in Fig. (b). The case where HEL ~ HCL (occurring for small barriers) is 


studied in more detail in section [TV] where we discuss how chemistry directly affect 
the different energy transferring mechanisms. We find that lowering the activation 
barrier on the chemical reactions leads to faster chemistry and a stronger coupling to 
the leading shock wave, and therefore, to a reduction in the pressure in this region. 
Whereas the propagation of the chemical and elastic waves achieve steady-state dur¬ 
ing the simulation, an accurate determination of the velocity plastic wave is difficult 
for the sizes of the system simulated as the plastic wave separates from the elastic 
precursor at late times, see Fig. (b). 

The effect of endothermicity was found to be minimal for all the cases considered 
in this study. Figure [TOj compares the effect of endothermicity on the Hugoniot and 
u^-Up curves for an inert and SWED materials with a curvature on the intramolecular 
potential corresponding to (AE/Vo)^dvi/ = 20 %. We see that, for cases where volume 
collapse is small and larger activation barriers, endothermicity can slightly alleviate 
the shock speed, and therefore, its resultant pressure. We note that while the an¬ 
thracene data matches the inert (AE/Vo)i;dm = 35 % case, in the (AE/Vo)i;dm = 20% 
case, the Hugoniot appears to match a reactive case with activation barriers between 
30 — 60 kcal/mol. 


IV. ROLE OF KINETICS AND LACK OF LOCAL EQUILIBRIUM IN THE 
NUCLEATION OF CHEMISTRY 

In the previous sections, we characterized the interplay between chemistry and 
the shock response for a family of SWED materials. We now focus on the molecular- 
level mechanisms that control the nucleation of volume-collapsing chemical reactions. 
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FIG. 10. Hugoniot and u^-Up curves for an Inert and SWED materials with a curvature 
on the intramolecular potential corresponding to /VQ)^dw = 20%. Different symbols 
correspond to different barriers. Empty symbols correspond no endothermicity {AH — 0 
kcal/mol) and filled symbols corresponds to AH — 20 kcal/mol 

We are interested, specifically, in understanding how the translational energy of the 
shockwave is transferred to the modes responsible for chemistry, i.e. the breathing 
mode, and how this process influences the initiation of chemistry. To do so, we study 
how the excitation energy from the shockwave is partitioned between the various 
meso, radial and internal modes and explore whether the coupling rates between 
them can influence chemical reactions. 

We focus on the response of a SWED material with AG = 30 kcal/mol and 
{AV/VQ)^dw = 35 %. Such a sample will be shocked at Up = 0.75 km/s as this 
corresponds to the HCL for this set of parameters. Figure pT] shows the time evolution 
of the temperatures associated with the three sets of DoFs for two thin slabs of 
material; one 20 nm away from the impact surface and another 60 nm away (inset 
in figure). The set of plots capture the effect of the coupling constants between the 
implicit DoFs and the particles c.m. [fmeso] and radial (or breathing) modes {urad)- 
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The meso temperature is further divided into normal (along the shock direction) and 
transverse modes. A dashed green line represent the percentage of reactive molecules 
as a function of time. Atomics snapshots of reacted particles and crystal structure 
are also shown at time 50 ps; fee is colored in green, bcc in blue, hep in white, and 
reacted particles in red. It is very clear that the rates associated with energy transfer 
affect chemistry. 

As the shock passes through the material, the molecular system will be away 
from local equilibrium with the modes that couple more strongly to the shock having 
higher temperatures. In all cases, the energy in the shockwave initially excites the 
particles c.m. and breathing modes as these modes are involved directly in the 
propagation of the compression wave. These modes achieve high temperatures in very 
short timescales and this excitation is transferred to the internal modes over longer 
timescales, which depend on the coupling rates, I'mess and I'rad- This is consistent 
with all atom MD simulation, see for example Refs. daEO]. 

The high temperature in the radial modes, together with the high pressure caused 
by the shockwave facilitates the volume-reducing chemical reactions. As the barrier 
for chemical reactions is overcome and the molecules relax to the metastable, low- 
volume state, the breathing modes are excited and cases with significant chemistry 
are marked by higher transient temperatures in the breathing modes. On the other 
hand, in regions with little chemistry, the mesoscale and radial temperatures achieve 
equilibrium faster. 

Interestingly, ChemDID predicts that chemical reactions can be initiated within 
short timescales, during this non-equilibrium stage. Thus, the timescales of energy 
exchange between the various sets of modes which determine the range of temper¬ 
ature experienced by the breathing modes play an important role in the chemical 
reactions. The nucleation of chemistry is also facilitated by high stresses and nu- 
cleation is observed predominantly in active slip planes where plastic deformation is 


27 


Temperature [K] Temperature [K] 


localized. This can be seen in the atomic snapshots shown in Fig. m 




(C) (D) 

FIG. 11. Atomic snapshots showing the chemical conversion (red particles) for various rates 
of the inter-molecular and intra-molecular ly^ad couplings. The plots below follow the 

time dependence of the temperature components for two thin slabs at the positions denoted 
by the arrows. Varying the rate of inter-molecular and intra-molecular couplings quenches 
the amount of chemistry (dashed line) at different rates as shown on the scale on the right 
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We now focus on the effect of energy coupling rates on the chemical response of the 


materials. To help in the discussion, Fig. [12] shows the pressure profiles for different 
coupling rates we have used in Fig{T^ For weak coupling rates to the internal DoFs, 
I'meso = I'rad = 0.01 ps“^ in Fig. 11 (a), the radial and mesoscopic temperatures 
remain in the non-equilibrium, high-temperature state for extended periods of time 
leading to significant chemical reactions. The propagation of a chemical front leads 
to a significant reduction in pressure behind the leading shockwave, see red line in 
Fig. As a result, this pressure reduction leads to a significant quenching of 

chemistry on the region behind the shock, except between active slip planes with 
high stress. 

Increasing the internal-to-intermolecular coupling {Pmeso = 0.1 ps“^) and removing 
the internal-to-radial coupling {Prad = 0.0 ps“^), leads to a reduction in chemical 
conversion throughout the sample. This occurs rather indirectly, as the internal 
temperature couples to the mesoscopic temperature via the DID equations and the 
latter couples to the radial DoFs via the Hamiltonian. The result is that the radial 
temperatures are lower compared to the case with weak coupling discussed above. 


see Fig. 11 (b). Now reversing the implicit couplings, i.e. Vrad — 1-0 ps“^ and 
^meso = 0.0 ps“^, leads to faster cooling on the radial temperature by the internal 
DoFs and a much slower chemical reaction front propagation. As expected, quenching 
the amount of chemistry near the impact leads to higher pressure in the region behind 
the leading shock, and as a result, to slightly more chemistry (4 %) in the region 
behind the leading shock, compared to the other two previous cases (<2%). 

Lastly, strong coupling between internal and meso and radial modes, r'meso = 
^rad = 0.1 ps“^, results in rapid equilibration of the radial and mesoscale tempera¬ 
tures with the internal temperature. We see equilibration with all the temperatures 
occurring within 30 ps throughout the whole sample, leading to a fast quenching 
effect on chemistry production. The pressure for this case, is the highest compared 
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FIG. 12. Pressure profiles at time = 50 ps after impact for various coupling values 

to all the other cases. As a result we see a large fraction of nucleation of chemical 
reaction behind the leading shock, but interestingly, these nucleations do not subse¬ 
quently grow into a propagating chemical wave. Our simulations clearly show that 
the ability of the chemical reactions to weaken shock waves is also affected by the 
details of the energy transfer between modes in the molecules. 


V. CONCLUSIONS 

In this paper we used ChemDID mesoscale simulations to explore the shock re¬ 
sponse of materials that can undergo endothermic, volume reducing chemical reac¬ 
tions. The simulations demonstrate that such stress-induced chemical reactions can 
be effective at dissipating the shockwave and such materials could be valuable in 
applications requiring protection from dynamical loads. 

We find that when a critical shock strength is reached, a chemical reaction wave 
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is induced behind the leading shock wave. The chemical reactions weaken the lead¬ 
ing shockwave. For a range piston velocities above this hugoniot chemical limit we 
observe a three-wave structure, with the elastic and plastic waves traveling at similar 
speeds and a the chemical wave trailing behind. Increasing the piston velocity results 
in an increase in the velocity of the chemical wave but does not affect the velocity 
leading shock nor the corresponding pressure. Under this regime, the material weak¬ 
ens the shock. This effect continues until the chemical wave moves at the same speed 
as the leading wave. In such an overdriven driven regime an increase in the shock 
velocity is translated directly to the leading shock. 

The results of the simulations together with an analysis of the shockwave struc¬ 
ture using conservation laws show that both the amount of volume collapse and the 
velocity of the chemical wave (governed by the activation energy associated with the 
chemical reactions) are the critical parameters for shockwave attenuation. The en- 
dothermicity of the chemical reactions plays a secondary role that is only discernible 
in the case of high barriers and modest volume collapse. Such information should be 
useful in the design and optimization of materials for shockwave energy dissipation. 

Interestingly, energy transfer rates between the various modes in the material 
also play a role in the nucleation and propagation of chemical reactions. Dynami¬ 
cal loading of materials leads to non-equilibrium states right behind the shockwave; 
with various DoFs experiencing different temperatures. Chemical reactions are in¬ 
fluenced by the temperature history of the corresponding modes and the coupling 
rates between them. Recent reactive MD simulations of solid explosives also indi¬ 
cate the possibility of chemical reactions away from local equilibrium following shock 
loading. ra The interplay between plasticity and chemistry will also affect resultant 
chemical reactions and shock attenuation in the material. Whereas our simulations 
presented focused on single crystalline materials, ChemDID could also be a useful 
tool to explore the effect of microstructure and defects on shock induced chemistry. 
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This would require large-scale simulations and we are working on an implementation 
of our mesoscale model in the LAMMPS parallel simulator that would enable such 
studies. 

For a given piston speed, a specific volume collapse provides a limited range 
where chemistry is able alleviate pressure buildup. Higher volume-collapse reduces 
the shock pressure the most, but also reduces the shock velocity corresponding to 
the aforementioned overdriven state. Avoiding this overdriven regime may be of 
interest for some applications. As a function of piston speeds, the larger activation 
energies can postpone the onset of chemistry to higher piston speeds before it reaches 
the overdriven state. Therefore, it is important to take both of these effects into 
consideration in order to design SWED materials that achieve the maximum amount 
of pressured dissipated over a desired range of velocities. 

In this paper we focused on the idealized case of sustained shocks and model 
materials. Additional complexity is introduced when considering finite pistons that 
lead to reflection and rarefaction waves that depend on the relative size of the target 
to impactor. In such cases non-steady states need to be considered, yet the general 
conclusions presented in this paper are useful for such cases. Finally, ChemDID can 
be parameterized to describe specific materials of interest. This would likely require 
more complex inter-molecular interaction potentials and chemical kinetics matching 
the real case. While developing such models is challenging, the technique would 
enable achieving time and length scales well beyond what is possible today with 
all-atom MD. 
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VI. APPENDIX 


A. Pressure quantification 


We write an expression for the pressure between a chemical region {c} and a 
shocked region {s} where none of the components are at rest. The transition across 
this region has to satisfy momentum conservation: 


Tc T Pc(^C ^c) Ps{,^s ^c) 

Let us expand this expression, we get 

Pc Ps Ps{,^s ^c) Pci,^c ^c) 

= Ps{xl + ul- 2XsUc) - Pc{xl +ul- 2XcUc) 

= PsXl - PcXl + Uc[ps{Uc - 2Xs) - PciUc - 2Xc)] 

PsXg PcXq T ^c[Ps(^c ^s) Pcip^c ^c) T PcXc 

= PgX^^ - PcXl + Uc[pcXc - PsXs] 

PcXcip^c Xc) T psXsiXs ^c) 


,where we have made use of mass conservation ( Eqn, 11 ) above. Substituting Eqn. 


12, twice, and rearranging, we obtain the final form 


Pc- Ps = (w - Xcfpc{— - 1) 

Ps 
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